##load data and packages

pacman::p_load(fixest, tidyverse, broom, data.table)

df <- readRDS("Data/Regimedata/Measures_merged_covariates.rds")
names(df)
summary(df$GDP_PC_Farris)


# normalize cagenmob and physical violence measure

df$v2cagenmob_nor <- (df$v2cagenmob - min(df$v2cagenmob, na.rm = T)) / (max(df$v2cagenmob, na.rm = T) - min(df$v2cagenmob, na.rm = T))
df$v2x_clphy*-1
df$v2x_clphy_nor <- (df$v2x_clphy - min(df$v2x_clphy, na.rm = T)) / (max(df$v2x_clphy, na.rm = T) - min(df$v2x_clphy, na.rm = T)) 

# select commnon sample in one df

df2 <- df %>% 
  select(RegType_lied_EA, RegType_RoW_EA,  Politytype_Anocracy, status_fh_PF,
         HTW_RegType_MP_Autocracy, RegType_magaloni_EA, AnckarRegtype_MP_Autocracy, year, Final_Code) %>% 
  na.omit() %>%
  select(year, Final_Code) %>%
  left_join(df, by = c("year", "Final_Code"))



## create dfs based on different typologies and run models with and without FE for the regime breakdown and the democratization outcomes

combinations <- tibble(Sample = c("RegType_lied_EA", "RegType_RoW_EA", "Politytype_Anocracy", "status_fh_PF", "AnckarRegtype_MP_Autocracy", "HTW_RegType_MP_Autocracy","RegType_magaloni_EA"),
                name = c("LIED", "RoW", "Polity", "FH", "CPR", "ARD", "AoW"))





results <- apply(combinations, 1, function(x) { 
Regime <- x[1]
name <- x[2]

list(

 feols(reg_breakdown ~ negative_growth | Final_Code + year, cluster = "Final_Code", data = df %>% filter (!!sym(Regime) == 1) ) %>%tidy()%>%
   mutate(Sample = "Full", Pred = "Recession", Measure = name),
 
 feols(reg_breakdown ~ log(GDP_PC_Farris) | Final_Code + year, cluster = "Final_Code", data = df %>% filter (!!sym(Regime) == 1) ) %>%tidy()%>%
   mutate(Sample = "Full", Pred = "log(GDP/PC)", Measure = name),
  
 feols(reg_breakdown ~ v2cagenmob| Final_Code + year, cluster = "Final_Code", data = df %>% filter ( !!sym(Regime)  == 1) ) %>%tidy()%>%
  mutate(Sample = "Full", Pred = "Mobilization", Measure = name),
  
 feols(reg_breakdown ~ v2x_clphy_nor | Final_Code + year, cluster = "Final_Code", data = df %>% filter ( !!sym(Regime)  == 1) ) %>%tidy()%>%
  mutate(Sample = "Full", Pred = "Repression", Measure = name),
  
 feols(reg_breakdown ~ negative_growth | Final_Code + year, cluster = "Final_Code", data = df2 %>% filter ( !!sym(Regime)  == 1) ) %>%tidy()%>%
   mutate(Sample = "Overlapping Sample", Pred = "Recession", Measure = name),
 
 feols(reg_breakdown ~ log(GDP_PC_Farris)  | Final_Code + year, cluster = "Final_Code", data = df2 %>% filter ( !!sym(Regime)  == 1) ) %>%tidy()%>%
   mutate(Sample = "Overlapping Sample", Pred = "log(GDP/PC)", Measure = name),
  
 feols(reg_breakdown ~ v2cagenmob| Final_Code + year, cluster = "Final_Code", data = df2 %>% filter ( !!sym(Regime)  == 1) ) %>%tidy() %>%
   mutate(Sample = "Overlapping Sample", Pred = "Mobilization", Measure = name),

  
  feols(reg_breakdown ~ v2x_clphy_nor | Final_Code + year, cluster = "Final_Code", data = df2 %>% filter ( !!sym(Regime)  == 1) )  %>%tidy() %>%
    mutate(Sample = "Overlapping Sample", Pred = "Repression", Measure = name)) %>% rbindlist()
  
})



coefplotframe <- results %>% 
  rbindlist %>% 
  as.data.frame %>% filter ( term != "(Intercept)" ) %>% mutate ( conf.low_95 = estimate - 1.96 *  std.error, conf.high_95 = estimate + 1.96 * std.error,
                                                                  conf.low_90 = estimate - 1.64 *  std.error, conf.high_90 = estimate + 1.64 * std.error)
coefplotframe$Measure <- fct_rev(factor(coefplotframe$Measure, levels=c("LIED", "RoW", "Miller", "Polity", "FH", "ARD", "AoW", "CPR")))



coefplotframe_bd <- coefplotframe %>% mutate(FE  = "Country & Year", Outcome = "Regime Breakdown")


results <- apply(combinations, 1, function(x) { 
  Regime <- x[1]
  name <- x[2]
  
  list(
    
    feols(reg_breakdown ~ negative_growth, cluster = "Final_Code", data = df %>% filter (!!sym(Regime) == 1) ) %>%tidy()%>%
      mutate(Sample = "Full", Pred = "Recession", Measure = name),
    
    feols(reg_breakdown ~ log(GDP_PC_Farris), cluster = "Final_Code", data = df %>% filter (!!sym(Regime) == 1) ) %>%tidy()%>%
      mutate(Sample = "Full", Pred = "log(GDP/PC)", Measure = name),
    
    feols(reg_breakdown ~ v2cagenmob, cluster = "Final_Code", data = df %>% filter ( !!sym(Regime)  == 1) ) %>%tidy()%>%
      mutate(Sample = "Full", Pred = "Mobilization", Measure = name),
    
    feols(reg_breakdown ~ v2x_clphy, cluster = "Final_Code", data = df %>% filter ( !!sym(Regime)  == 1) ) %>%tidy()%>%
      mutate(Sample = "Full", Pred = "Repression", Measure = name),
    
    feols(reg_breakdown ~ negative_growth, cluster = "Final_Code", data = df2 %>% filter ( !!sym(Regime)  == 1) ) %>%tidy()%>%
      mutate(Sample = "Overlapping Sample", Pred = "Recession", Measure = name),
    
    feols(reg_breakdown ~ log(GDP_PC_Farris), cluster = "Final_Code", data = df2 %>% filter ( !!sym(Regime)  == 1) ) %>%tidy()%>%
      mutate(Sample = "Overlapping Sample", Pred = "log(GDP/PC)", Measure = name),
    
    
    feols(reg_breakdown ~ v2cagenmob, cluster = "Final_Code", data = df2 %>% filter ( !!sym(Regime)  == 1) ) %>%tidy() %>%
      mutate(Sample = "Overlapping Sample", Pred = "Mobilization", Measure = name),
    
    
    feols(reg_breakdown ~ v2x_clphy, cluster = "Final_Code", data = df2 %>% filter ( !!sym(Regime)  == 1) )  %>%tidy() %>%
      mutate(Sample = "Overlapping Sample", Pred = "Repression", Measure = name)) %>% rbindlist()
  
})




coefplotframe <- results %>% 
  rbindlist %>% 
  as.data.frame %>% filter ( term != "(Intercept)" ) %>% mutate ( conf.low_95 = estimate - 1.96 *  std.error, conf.high_95 = estimate + 1.96 * std.error,
                                                                  conf.low_90 = estimate - 1.64 *  std.error, conf.high_90 = estimate + 1.64 * std.error)
coefplotframe$Measure <- fct_rev(factor(coefplotframe$Measure, levels=c("LIED", "RoW", "Miller", "Polity", "FH", "ARD", "AoW", "CPR")))

coefplotframe_bd2 <- coefplotframe %>% mutate(FE  = "No FE", Outcome = "Regime Breakdown")







combinations <- tibble(Sample = c("RegType_lied_EA", "RegType_RoW_EA", "Politytype_Anocracy", "status_fh_PF", "AnckarRegtype_MP_Autocracy", "HTW_RegType_MP_Autocracy","RegType_magaloni_EA"),
                       name = c("LIED", "RoW", "Polity", "FH", "CPR", "ARD", "AoW"),
                       Depvar =  c("F_RegType_lied_DEM", "F_RegType_RoW_DEM", "F_Politytype_Democracy", "F_status_fh_F", "F_AnckarRegtype_Democracy", "F_HTW_RegType_Democracy", "F_RegType_magaloni_DEM"))


names(df)




results2 <- apply(combinations, 1, function(x) { 
  Regime <- x[1]
  name <- x[2]
  depvar <- x[3]
  
  list(
    
    feols(as.formula(paste(depvar, "~", "log(GDP_PC_Farris)", "| Final_Code + year")), cluster = "Final_Code", data = df %>% filter (!!sym(Regime) == 1) ) %>%tidy()%>%
      mutate(Sample = "Full", Pred = "log(GDP/PC)", Measure = name),
    
    feols(as.formula(paste(depvar, "~", "negative_growth", "| Final_Code + year")), cluster = "Final_Code", data = df %>% filter (!!sym(Regime) == 1) ) %>%tidy()%>%
      mutate(Sample = "Full", Pred = "Recession", Measure = name),
    
    feols(as.formula(paste(depvar, "~", "v2cagenmob", "| Final_Code + year")), cluster = "Final_Code", data = df %>% filter ( !!sym(Regime)  == 1) ) %>%tidy()%>%
      mutate(Sample = "Full", Pred = "Mobilization", Measure = name),
    
    feols(as.formula(paste(depvar, "~", "v2x_clphy", "| Final_Code + year")), cluster = "Final_Code", data = df %>% filter ( !!sym(Regime)  == 1) ) %>%tidy()%>%
      mutate(Sample = "Full", Pred = "Repression", Measure = name),
    
    feols(as.formula(paste(depvar, "~", "reg_breakdown", "| Final_Code + year")), cluster = "Final_Code", data = df %>% filter ( !!sym(Regime)  == 1) ) %>%tidy()%>%
      mutate(Sample = "Full", Pred = "Reg. Breakdown", Measure = name),
    
    
    feols(as.formula(paste(depvar, "~", "log(GDP_PC_Farris)", "| Final_Code + year")), cluster = "Final_Code", data = df2 %>% filter ( !!sym(Regime)  == 1) ) %>%tidy()%>%
      mutate(Sample = "Overlapping Sample", Pred = "log(GDP/PC)", Measure = name),
    
    feols(as.formula(paste(depvar, "~", "negative_growth", "| Final_Code + year")), cluster = "Final_Code", data = df %>% filter (!!sym(Regime) == 1) ) %>%tidy()%>%
      mutate(Sample = "Overlapping Sample", Pred = "Recession", Measure = name),
    
    
    feols(as.formula(paste(depvar, "~", "v2cagenmob", "| Final_Code + year")), cluster = "Final_Code", data = df2 %>% filter ( !!sym(Regime)  == 1) ) %>%tidy() %>%
      mutate(Sample = "Overlapping Sample", Pred = "Mobilization", Measure = name),
    
    feols(as.formula(paste(depvar, "~", "v2x_clphy", "| Final_Code + year")), cluster = "Final_Code", data = df2 %>% filter ( !!sym(Regime)  == 1) )  %>%tidy() %>%
      mutate(Sample = "Overlapping Sample", Pred = "Repression", Measure = name),
    
    feols(as.formula(paste(depvar, "~", "reg_breakdown", "| Final_Code + year")), cluster = "Final_Code", data = df2 %>% filter ( !!sym(Regime)  == 1) )  %>%tidy() %>%
      mutate(Sample = "Overlapping Sample", Pred = "Reg. Breakdown", Measure = name)) %>% rbindlist()
  
})




coefplotframe <- results2 %>% 
  rbindlist %>% 
  as.data.frame %>% filter ( term != "(Intercept)" ) %>% mutate ( conf.low_95 = estimate - 1.96 *  std.error, conf.high_95 = estimate + 1.96 * std.error,
                                                                  conf.low_90 = estimate - 1.64 *  std.error, conf.high_90 = estimate + 1.64 * std.error)
coefplotframe$Measure <- fct_rev(factor(coefplotframe$Measure, levels=c("LIED", "RoW", "Miller", "Polity", "FH", "ARD", "AoW", "CPR")))


coefplotframe_de <- coefplotframe %>% mutate(FE  = "Country & Year", Outcome = "Democratization")





results3 <- apply(combinations, 1, function(x) { 
  Regime <- x[1]
  name <- x[2]
  depvar <- x[3]
  
  list(
    
    feols(as.formula(paste(depvar, "~", "log(GDP_PC_Farris)")), cluster = "Final_Code", data = df %>% filter (!!sym(Regime) == 1) ) %>%tidy()%>%
      mutate(Sample = "Full", Pred = "log(GDP/PC)", Measure = name),
    
    feols(as.formula(paste(depvar, "~", "negative_growth")), cluster = "Final_Code", data = df %>% filter (!!sym(Regime) == 1) ) %>%tidy()%>%
      mutate(Sample = "Full", Pred = "Recession", Measure = name),
    
    feols(as.formula(paste(depvar, "~", "v2cagenmob")), cluster = "Final_Code", data = df %>% filter ( !!sym(Regime)  == 1) ) %>%tidy()%>%
      mutate(Sample = "Full", Pred = "Mobilization", Measure = name),
    
    feols(as.formula(paste(depvar, "~", "v2x_clphy")), cluster = "Final_Code", data = df %>% filter ( !!sym(Regime)  == 1) ) %>%tidy()%>%
      mutate(Sample = "Full", Pred = "Repression", Measure = name),
    
    feols(as.formula(paste(depvar, "~", "reg_breakdown")), cluster = "Final_Code", data = df %>% filter ( !!sym(Regime)  == 1) ) %>%tidy()%>%
      mutate(Sample = "Full", Pred = "Reg. Breakdown", Measure = name),
    
    
    feols(as.formula(paste(depvar, "~", "log(GDP_PC_Farris)")), cluster = "Final_Code", data = df2 %>% filter ( !!sym(Regime)  == 1) ) %>%tidy()%>%
      mutate(Sample = "Overlapping Sample", Pred = "log(GDP/PC)", Measure = name),
    
    feols(as.formula(paste(depvar, "~", "negative_growth")), cluster = "Final_Code", data = df2 %>% filter ( !!sym(Regime)  == 1) ) %>%tidy()%>%
      mutate(Sample = "Overlapping Sample", Pred = "Recession", Measure = name),
    
    feols(as.formula(paste(depvar, "~", "v2cagenmob")), cluster = "Final_Code", data = df2 %>% filter ( !!sym(Regime)  == 1) ) %>%tidy() %>%
      mutate(Sample = "Overlapping Sample", Pred = "Mobilization", Measure = name),
    
    feols(as.formula(paste(depvar, "~", "v2x_clphy")), cluster = "Final_Code", data = df2 %>% filter ( !!sym(Regime)  == 1) )  %>%tidy() %>%
      mutate(Sample = "Overlapping Sample", Pred = "Repression", Measure = name),
    
    feols(as.formula(paste(depvar, "~", "reg_breakdown")), cluster = "Final_Code", data = df2 %>% filter ( !!sym(Regime)  == 1) )  %>%tidy() %>%
      mutate(Sample = "Overlapping Sample", Pred = "Reg. Breakdown", Measure = name)) %>% rbindlist()
  
})

names(df)


coefplotframe <- results3 %>% 
  rbindlist %>% 
  as.data.frame %>% filter ( term != "(Intercept)" ) %>% mutate ( conf.low_95 = estimate - 1.96 *  std.error, conf.high_95 = estimate + 1.96 * std.error,
                                                                  conf.low_90 = estimate - 1.64 *  std.error, conf.high_90 = estimate + 1.64 * std.error)
coefplotframe$Measure <- fct_rev(factor(coefplotframe$Measure, levels=c("LIED", "RoW", "Miller", "Polity", "FH", "ARD", "AoW", "CPR")))


coefplotframe_de2 <- coefplotframe %>% mutate(FE  = "No FE", Outcome = "Democratization")



coefplotframe <- rbind(coefplotframe_bd, coefplotframe_bd2,
                       coefplotframe_de, coefplotframe_de2)


coefplotframe <- coefplotframe %>% 
  filter(Measure != "Miller")

coefplotframe$Measure <- as.factor(as.character(coefplotframe$Measure))


levels(coefplotframe$Measure) <- c("Authoritarian Regimes Data Set", "Autocracies of the World", "Classifying Political Regimes", "Freedom House", "Lexical Index of Electoral Democracy", "Polity", "Regimes of the World")

coefplotframe$Measure <- factor(coefplotframe$Measure, levels=c("Lexical Index of Electoral Democracy", "Regimes of the World", "Polity", "Freedom House", "Authoritarian Regimes Data Set", "Autocracies of the World", "Classifying Political Regimes"))



### Plots results

uniqe(coefplotframe$Pred)


coefplotframe$Pred <- as.factor(coefplotframe$Pred)
levels(coefplotframe$Pred)
coefplotframe$Pred <- factor(coefplotframe$Pred  , levels=c("Mobilization", "Recession", "Reg. Breakdown", "Repression", "log(GDP/PC)") )

### Figure G1

coefplotframe %>% filter(FE == "Country & Year" ) %>% 
  ggplot(aes(x=Sample, y = estimate, shape = fct_rev(Measure))) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_linerange(aes(x = Sample, 
                     ymin = conf.low_95,
                     ymax = conf.high_95), position = position_dodge(width = .75), alpha = .7, color = "grey40") +
  geom_linerange(aes(x = Sample, 
                     ymin = conf.low_90,
                     ymax = conf.high_90), position = position_dodge(width = .75), linewidth = 1, alpha = .7, color = "grey40")   +
  geom_point(aes(x = Sample, 
                 y = estimate), size = 4, position = position_dodge(width = .75), color = "black") +
  ggtitle("") + ylab ("") + xlab("") +  guides(color = guide_legend(reverse=TRUE, ncol = 2, title = ""), shape = guide_legend(reverse=TRUE, ncol = 2, title = "")) + 
  coord_flip()  + scale_y_continuous(breaks = 0) +
  theme_classic(base_size = 20) + facet_grid(vars(fct_rev(Outcome)), vars(fct_rev(Pred)),  scales = "free_x") +
  labs(colour="Measure", shape = "Measure")  + scale_color_grey() +
  scale_shape_manual(values = c(55:49)) + theme(legend.position="bottom")


### Figure H2

coefplotframe %>% filter(FE == "No FE" ) %>%   ggplot(aes(x=Sample, y = estimate, shape = fct_rev(Measure))) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_linerange(aes(x = Sample, 
                     ymin = conf.low_95,
                     ymax = conf.high_95), position = position_dodge(width = .75), alpha = .7, color = "grey40") +
  geom_linerange(aes(x = Sample, 
                     ymin = conf.low_90,
                     ymax = conf.high_90), position = position_dodge(width = .75), linewidth = 1, alpha = .7, color = "grey40")   +
  geom_point(aes(x = Sample, 
                 y = estimate), size = 4, position = position_dodge(width = .75), color = "black") +
  ggtitle("") + ylab ("") + xlab("") +   guides(color = guide_legend(reverse=TRUE, ncol = 2, title = ""), shape = guide_legend(reverse=TRUE, ncol = 2, title = "")) + 
  coord_flip()  + scale_y_continuous(breaks = 0) +
  theme_classic(base_size = 20) + facet_grid(vars(fct_rev(Outcome)), vars(fct_rev(Pred)),  scales = "free_x") +
  labs(colour="Measure", shape = "Measure")  + scale_color_grey() +
  scale_shape_manual(values = c(55:49))+ theme(legend.position="bottom")





## Figure I2

coefplotframe %>% filter(FE == "Country & Year" ) %>% 
  ggplot(aes(x=Sample, y = estimate, shape = fct_rev(Measure))) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_linerange(aes(x = Sample, 
                     ymin = conf.low_95,
                     ymax = conf.high_95), position = position_dodge(width = .75), alpha = .7, color = "grey40") +
  geom_linerange(aes(x = Sample, 
                     ymin = conf.low_90,
                     ymax = conf.high_90), position = position_dodge(width = .75), linewidth = 1, alpha = .7, color = "grey40")   +
  geom_point(aes(x = Sample, 
                 y = estimate), size = 4, position = position_dodge(width = .75), color = "black") +
  ggtitle("") + ylab ("") + xlab("") +  guides(color = guide_legend(reverse=TRUE, ncol = 2, title = ""), shape = guide_legend(reverse=TRUE, ncol = 2, title = "")) + 
  coord_flip()  + scale_y_continuous(breaks = c(-.5,0,.5)) +
  theme_classic(base_size = 20) + facet_grid(vars(fct_rev(Outcome)), vars(fct_rev(Pred))) +
  labs(colour="Measure", shape = "Measure")  + scale_color_grey() +
  scale_shape_manual(values = c(55:49)) + theme(legend.position="bottom")


